#####################################################################
# Table 3
#####################################################################

rm(list=ls())

library(lfe)
library(Hmisc)
library(ggplot2)
library(reshape2)
library(xtable)
library(stargazer)
library(lmtest)
library(sandwich)
library(stargazer)
library(msm)

# load function that does clustered SEs
vcovCluster <- function(
  model,
  cluster
)
{
  require(sandwich)
  require(lmtest)
  if(nrow(model.matrix(model))!=length(cluster)){
    stop("check your data: cluster variable has different N than model")
  }
  M <- length(unique(cluster))
  N <- length(cluster)           
  K <- model$rank   
  if(M<40){
    warning("Fewer than 40 clusters, variances may be unreliable")
  }
  dfc <- (M/(M - 1)) * ((N - 1)/(N - K))
  uj  <- apply(estfun(model), 2, function(x) tapply(x, cluster, sum));
  rcse.cov <- dfc * sandwich(model, meat = crossprod(uj)/N)
  return(rcse.cov)
}

############################
# Read and prepare data 
############################

# prepare and load data from LAPOP website 

# drop missing values
d0 = d
d = na.omit(d)

##########################
# Adjusted results
##########################

mod1 = felm(d$crossline ~ d$victimization + d$education + d$female + d$age  + as.factor(d$ethnicity) | as.factor(d$country) | 0 | d$municipality)
summary(mod1)

mod2 = felm(d$supportdemocracy ~ d$victimization + d$education + d$female + d$age + as.factor(d$ethnicity)  | as.factor(d$country) | 0 | d$municipality)
summary(mod2)

##############
# Tables
##############

# Table 1
stargazer(mod1,mod2,
          type = "text",
          title="Regression results",
          align=TRUE, 
          omit.stat=c("LL","ser","f","rsq","adj.rsq"), 
          no.space=TRUE,  
          multicolumn = TRUE,
          table.placement = "H",
          #single.row = TRUE,
          #column.separate = c(2,2),
          covariate.labels=c("Crime victimization"),
          #dep.var.caption  = "Strong-handed policies",
          column.labels = c("Strong-handed policies","Support for democracy"),
          dep.var.labels.include = FALSE,
          #star.cutoffs = c(0.05,0.01,0.001),
          omit = c("education","female","age","ethnicity","province","year","Constant"),
          add.lines = list(c("Placebo covariates","Yes","Yes"),
                           c("Country fixed effects","Yes","Yes"),
                           c("Countries","18","18")))
